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Diffusion tensor imaging (DTI) plays a key role in analyzing the 

physical structures of biological tissues, particularly in reconstructing 

fiber tracts of the human brain in vivo. On the one hand, eigenvalues 

of diffusion tensors (DTs) estimated from diffusion weighted imag- 

Mh , ing (DWI) data usually contain systematic bias, which subsequently 

biases the diffusivity measurements popularly adopted in fiber track- 
ing algorithms. On the other hand, correctly accounting for the spa- 
C^ , tial information is important in the construction of these diffusiv- 

f/3 ' ity measurements since the fiber tracts are typically spatially struc- 

tured. This paper aims to establish test-based approaches to identify 
anisotropic water diffusion areas in the human brain. These areas 
in turn indicate the areas passed by fiber tracts. Our proposed test 
^^ , statistic not only takes into account the bias components in eigenvalue 

^^ ■ estimates, but also incorporates the spatial information of neighbor- 

r*--, ' ing voxels. Under mild regularity conditions, we demonstrate that 

^^ , the proposed test statistic asymptotically follows a x^ distribution 

under the null hypothesis. Simulation and real DTI data examples 
f— ^ \ are provided to illustrate the efficacy of our proposed methods. 
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1. Introduction. Diffusion tensor imaging (DTI) has been widely used 
by neuroscientists to reconstruct the pathways of white matter fibers in hu- 
man brain in vivo. DTI data are usually estimated from diffusion weighted 
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Fig. 1. Ellipsoid representation of a DT. 

imaging (DWI) data acquired in magnetic resonance experiments by some 
statistical model (see Section 2.1). A set of DTI data is typically com- 
posed of diffusion tensors (DTs), each contained in a corresponding voxel. 
Here, a voxel stands for a volume element in a 3D imaging space. Each 
DT, denoted by D, can be represented by a 3 x 3 symmetric positive def- 
inite matrix, which together with its decomposed eigenvalue-eigenvector 
pairs {(A(fc),V(fc)) : A(3) > A(2) > ^{i),k = 1,2,3} geometrically characterizes 
the degree and orientation of the water diffusion in that particular voxel. 
More specifically, the eigenvectors and the square root of eigenvalues of D, 
respectively, correspond to the orientations and the lengths of axes in an 
ellipsoid representation (see Figure 1). The distance between the center and 
any point on the surface of the ellipsoid measures the rate of the water 
diffusion along that particular orientation in the voxel. 

One of the main themes in DTI research is to identify the anisotropic 
water diffusion brain areas, which facilitates the downstream fiber tracking 
process. There are two general strategies aiming to address this problem. 
The first is to construct scalar measurements. The anisotropic water diffu- 
sion areas are then identified based on thresholding these scalar measure- 
ments. Reviews of this type of method can be found in Moseley et al. (1990), 
Douek et al. (1991), van Gelderen et al. (1994), Basser and Pierpaoli (1996), 
Johansen-Berg and Behrens (2009), among others. Thresholding the frac- 
tional anisotropy (FA) and the relative anisotropy (RA) [see Section 2.3; 
Basser and Pierpaoli (1996)] has gained popularity and has been widely 
adopted by neuroscientists in the past decade. Nonetheless, FA and RA are 
essentially defined as functions of the eigenvalues of the DT estimate in every 
brain voxel. These eigenvalues typically carry systematic bias [see Section 
3; Pierpaoli and Basser (1996), Zhu et al. (2007), Jones (2003), Lazar and 
Alexander (2003)], the magnitudes of which are sensitive to the distribution 
of the noise carried in raw DWI data, and therefore significantly affect the 
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effectiveness and validity of FA (RA) based methods in practice. The sec- 
ond is to classify the morphologies of DTs via test-based approaches [Zhu 
et al. (2006), Zhu et al. (2007)], which not only quantify the degree of wa- 
ter diffusivity in each voxel, but also lend theoretical supports to statistical 
inference. However, to the best of our knowledge, all existing approaches on 
this respect are single-voxel based. The validity and performance of these 
methods rely essentially on the technical requirement that the number of 
diffusion gradients in the DWI experiment is large and ultimately diverg- 
ing to infinity. Moreover, the DTI data are typically spatially structured. 
Ignoring the spatial information may diminish the effectiveness of the meth- 
ods. 

In this paper, we develop new test-based approaches to identify the aniso- 
tropic water diffusion brain areas, which are usually associated with areas 
passed by fiber tracts. To this end, for each voxel, we examine the testing 
problem "a// three eigenvalues are equivalent" against "ai least two eigen- 
values are different." The former corresponds to isotropically diffused DTs, 
whereas the latter includes anisotropically diffused DTs with morphologies 
of prolate, oblate and nondegenerate. Our proposed test statistic accom- 
modates the spatial information of the imaging space by taking into ac- 
count eigenvalues in neighboring voxels. Under mild regularity conditions, 
we demonstrate that our proposed test statistic asymptotically follows a x^ 
distribution. Therefore, the performance of our methods in the identification 
of fiber areas is not affected by the bias components carried in the eigen- 
value estimates. In theory, one of the main technical requirements is the 
divergence of the number of neighboring voxels involved in the construction 
of the statistic. This differs from the divergence of the number of diffusion 
gradients typically assumed by test-based approaches. Therefore, our meth- 
ods shed light on alternative ways of improving the identification accuracy 
of anisotropic water diffusion areas. Furthermore, an adaptive procedure to 
select varied neighborhoods is proposed to solidify the performance of our 
proposed approaches when the acquired imaging data have limited resolu- 
tion. Simulation studies and real data examples are provided to illustrate 
the efficacy of our proposed methods. 

The rest of the paper is organized as follows. Section 2 introduces the 
background related to our study. Section 3 establishes a statistical model 
based on eigenvalues in the selected neighborhood of a single voxel. Sec- 
tion 4 describes the procedure of constructing our proposed test statistic, 
and an adaptive method for selecting neighboring voxels. Section 5 explores 
the theoretical properties of our proposed test statistic. Section 6 presents 
simulation results. There, our methods are compared with FA-threshold and 
Smooth- FA-threshold approaches. Section 7 applies all approaches on real 
brain DTI data. Section 8 discusses our findings in this paper. Technical 
conditions and proofs are given in a supplemental document. 
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2. Background. We begin with a brief introduction of DWI data, DTI 
data and existing statistical models for estimating DTI data from DWI data. 
Then, we summarize the associations among fiber tracts, tissue types, water 
diffusivity and DT types. After that, we overview the quantitative scalars, 
FA and RA, popularly used in fiber tracking algorithms. 

2.1. From DWI to DTI. In this section we first give a brief introduction 
of the structures of DWI and DTI data, where the former are acquired 
from the diffusion weighted magnetic resonance experiment, while the latter 
are estimated from the former based on some statistical model. Then, we 
summarize the existing statistical models for estimating DTI data from DWI 
data. 

We assume that the DWI data over the brain of a given subject contain 
N voxels, each of which consists of diffusion-weighted measurements. De- 
note by (po, b, {(</'«, gj )}[=!, the acquired diffusion-weighted measurements 
at a given voxel over the brain in a DWI experiment. Here, the ith diffu- 
sion gradient gj = {gi,i,gi,2,9i,zf ■, with gfg, = 1, is chosen by the exper- 
imenter before the DWI experiment starts, and serves as a scanning di- 
rection in the experiment [Hasan, Parker and Alexander (2001)]; h is the 
6- factor, whose value is determined by a function of parameter settings in 
the DWI experiment [Stejskal and Tanner (1965), Moseley et al. (1990), 
Anderson (2001)]; both gj and h usually adopt the same values over all 
voxels in a DWI experiment; (f)i denotes the diffusion attenuated signal, 
acquired on the ith diffusion gradient gj at h] (pQ is the reference signal 
obtained at 6 = 0; {(/)j}[^^ and ^o compose the responses of the DWI exper- 
iment for each voxel; r is the number of acquired attenuated signals for each 
voxel. 

Accordingly, a single-voxel of the DTI data contains a 3 x 3 symmetric, 
positive definite DT matrix, 



D 



which carries the intrinsic information of water diffusion in that particu- 
lar voxel. The elements of D can be reorganized as a 6 x 1 vector d = 

iDi,l,D2,2,D3,3,Di^2,Di-s,D2,3f. 

The connections between DWI and DTI data are first investigated by 
the seminal work of Basser, Mattiello and LeBihan (1994), in which the 
following multivariate linear and nonlinear regression models are proposed. 
For a single-voxel, 

(2.1) multivariate linear model: log((;^j) = log{(f)o) — 6xj d + Sj, 

(2.2) multivariate nonlinear model: 0j = </)oexp(— 6xj d) + rji, 
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Fig. 2. An illustrative example of fiber tracts. 

where z = l,...,r, Xj = {gl^, gf^, gl.^,2gi^igi^2,'^gi,igi,3,'^9i,29i,3V , ei and r/j 
are random errors, d is then estimated from either model by regression 
techniques. 

A number of alternative statistical models as well as fitting procedures, 
besides Basser, Mattiello and LeBihan (1994), have been proposed to obtain 
sophisticated DT estimates from DWI data concerning various aspects, such 
as robustness, bias, non-Gaussian errors, spatial smoothness, model validity, 
etc. Examples include Mangin et al. (2002), Chang, Jones and Pierpaoli 
(2005), Salvador et al. (2005), Heim et al. (2007), Zhu et al. (2007), Tabelow 
et al. (2008) and many others. 



2.2. Fiber tracts, tissue types, water diffusivity and DT types. The as- 
sociations among fiber tracts, tissue types, water diffusivity and DT types 
are summarized as follows. For voxels located in fiber tracts, that is, white 
matter brain areas, water tends to present higher diffusivity along the domi- 
nant orientation of fibers than that in other orientations. DTs in these voxels 
are anisotropic, characterized by the heterogeneity in lengths of axes in the 
corresponding ellipsoid representations. In contrast, for voxels in brain areas 
without fiber tracts, that is, grey matter areas, DTs are isotropic. In these 
areas, the diffusivity of water in all orientations is roughly the same. There- 
fore, the corresponding represented ellipsoids are spherically shaped. The 
morphology of an anisotropic DT can usually be classified into one of the 
three categories, namely, prolate, oblate and nondegenerate, respectively, 
corresponding to brain voxels located in uniquely orientated fiber tracts 
{Area 1 in Figure 2), crossed fiber tracts with similar intensities on two 
or more different orientations {Area 2 in Figure 2, characterized by similar 
fiber denseness for the red and blue bundles in their intersected parts) and 
crossed fiber tracts with distinct intensities on different orientations {Area 3 
in Figure 2, characterized by the scenario that the denseness of blue tracts 
is higher than that of green tracts in their intersected parts) . 
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A vast number of tractography algorithms have been proposed in order to 
reconstruct fiber tracts in the human brain based on DTI data. A hst of ex- 
amples of these algorithms can be found in Conturo et al. (1999), Gossl et al. 
(2002), Xu et al. (2002), Behrens et al. (2007), O'Donneh and Westin (2007) 
and many others. We observe that most of these approaches are founded on 
the derived voxel-wise scalar quantities, such as FA and RA, which summa- 
rize/extract microstructural information of water diffusion carried by DWI 
or DTI data. 

2.3. Quantitative scalars: FA and RA. The fractional anisotropy (FA) 
and relative anisotropy (RA) [Basser and Pierpaoli (1996)] are quantitative 
scalars widely used by neuroscientists to measure the water diffusivity in 
brain tissues and construct algorithms for tracking fibers, since they are 
computationally simple and invariant in the choice of the laboratory coor- 
dinate system and diffusion gradients. For a given voxel, FA and RA are 
defined as 



FA 



3ELi{Aw-A(.)P ^^_ 3 ^Et=i{\k)-Xi.)V 



^3 /\,., _ \,.\2 



y ELiA?,) ' V2 E 
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where A(.) = {A(i) -I- A(2) + A(3)}/3, with A(3) > A(2) ^ ^{i) being the ordered 
eigenvalues of the DT estimate. 

Under the ideal but unrealistic assumption that the estimated DT is 
noise free, RA = and FA = in isotropic voxels [i.e., A(i) = A(2) = A(3)], 
whereas RA = -v/3 and FA = 1 in purely anisotropic voxels [i.e., A(3) ^> 
A(2) = A(i)]. However, the noise carried by the DWI data contaminates the 
DT estimates, and subsequently introduces systematic bias into the de- 
rived eigenvalues [Pierpaoli and Basser (1996)]. Although these bias compo- 
nents have been investigated by numerical evaluations [Jones (2003), Lazar 
and Alexander (2003)] as well as in theory [Zhu et al. (2007)], the mag- 
nitudes are sensitive to the distribution of the noise in DWI experiments. 
Consequently, these bias components introduce uncertainty into the con- 
structed FA and RA. In practice, brain areas with small but nonvanish- 
ing FA (RA) usually correspond to grey matter areas, whereas those with 
large FA (RA) are typically areas passed by fiber tracts. Some tractogra- 
phy algorithms are based on FA or RA with thresholds (e.g., the tractog- 
raphy algorithm integrated in MedlNRIA, a publicly available software at 
http : //www-sop . inria.fr/asclepios/software/MedINRIA). The thresh- 
olds, however, are usually manually chosen by investigators based on their 
historical knowledge of the DTI data and the structure of the human brain. 

Figure 3 illustrates how the fiber tracking results are affected by distinct 
experiential thresholds of the FA, where the detailed data information is 
given in Section 7. We use a region of interest (ROI) with 30 x 30 voxels. 
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Fig. 3. Fiber tracking results by distinct FA thresholds. From left to right: FA map with 
a ROI highlighted; constructed fibers passing through the ROI by FA thresholds 0.3, 0.45 
and 0.6, respectively. 



which is highhghted by a red rectangle shown in the leftmost panel of Fig- 
ure 3. The tracked fibers passing through the ROI using distinct FA thresh- 
olding criteria, namely, FA > c with c = 0.3, 0.45 and 0.6, are displayed in 
the remaining three panels of Figure 3. All sub-figures in Figure 3 are con- 
structed by MedlNRIA, where different colors stand for different principal 
orientations of the corresponding DTs. We observe that the fiber tracking 
results are sensitive to the FA threshold, the choice of which, owing to the 
uncertainty of the bias magnitudes in eigenvalue estimates, is not well sup- 
ported in theory. In other words, there is lack of a criterion to choose the 
threshold for FA, and justify the goodness of the tracking results. 

A more sophisticated approach based on FA is due to Zhu et al. (2006), 
in which the asymptotic distribution of a test statistic established on FA is 
investigated. Therefore, it is capable of suggesting a threshold to be used. 
This approach, however, is not applicable in our study, since the theory 
needs the assumption that the number of gradients for every brain voxel is 
diverging to infinity, that is, r — )■ oo. In contrast, in our study r = 12. 

Another limitation of FA and RA with threshold approaches is that the 
DTI data are typically spatially structured. FA and RA, however, are defined 
as functions of the eigenvalues of the DT in every single voxel. These ap- 
proaches identify the existence of fibers in each voxel while ignoring the spa- 
tial information, and therefore may diminish the effectiveness in the down- 
stream fiber tracking algorithms. 

We observe that there exist approaches which accommodate the spatial 
information in the stage of DT estimation, such as Heim et al. (2007) and 
Tabelow et al. (2008). Intuitively, by incorporating the spatial information, 
these methods improve the DT estimation, and therefore the accuracy of FA 
and RA. Hereafter, we refer to the FA based on DT estimated by the ap- 
proach in Tabelow et al. (2008) and Polzehl and Tabelow (2009) (implement 
in R package dti) as Smooth-FA. In our numerical studies. Smooth- FA with 
threshold, namely, Smooth-FA-threshold, is employed as an approach to iden- 
tify anisotropic brain voxels, and is compared with our proposed methods. 
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3. Statistical model on local eigenvalues. We consider the DTI data of 
a given subject with N voxels. In each voxel v € {1,. . . , N}, denote by D(u) 
(without confusion, we concisely denote D) the DT estimated from a statis- 
tical model in Section 2.1. Let A(3) > A(2) > ^(i) be the ordered eigenvalues 
of D. In practice, these eigenvalues are adopted to estimate the ordered true 
eigenvalues, denoted by ALx > \%-, > A^w 

It has been demonstrated both in numerical studies and in theory that 
£J{A(3)} > A^gN and £'{A(i)} < Xl^y Such a bias is actually caused by the 
sorting procedure in the decomposition of D. In other words, we can pos- 
tulate that one of the eigenvalues of D, that is, Afc G {^(k) '-k = 1,2,3}, is 
associated with AJ^n (but we don't know which one it is). Here, "associated" 
means there exists a random error e^ with E[ek) = such that 

(3.1) Afc = A^,)+efc. 

However, A/'^n is estimated by A(fc) instead of A^. Therefore, ii^{A(3)} > 
-^(As) = ATg-j. Similar arguments lead to £^{A(i)} < A^-,. The magnitudes 
of the bias components are sensitive to the distribution of the noise in DWI 
experiments, introducing uncertainty into the quantitative scalars, such as 
FA and RA. 

In our approach, instead of using the eigenvalues at a single voxel v only, 
we choose a set of n neighboring voxels located adjacent to v. Denote by 
{Aj (;.)(?;) : A; = 1, 2, 3}"^]^ [i.e., Xj^(k)\ the eigenvalue estimates in the selected 
neighboring voxels, whose permuted version associated with the set of true 
eigenvalues is denoted by {\j^k{v)'-k = 1,2,3}'^^^ (i.e., Xj^k)- Here, j = 1 
corresponds to voxel v, that is, Ai^fc = Afc. 

Similar in spirit to the random complete block design, we model E{\j^k), 
the expected eigenvalue of the jth neighboring voxel of v, as the addition of 
two components, namely, the true eigenvalue A/^-, for voxel v and the differ- 
ence of eigenvalues between voxel v and its jth neighboring voxel, denoted 
by (3j. That is, E{\j^k) = Kj.\ + /3j, which leads to the additive model, 

(3.2) Xj^k = A(;.) -I- /3j -I- Ej^k, 

where the eigenvalues {Aj!'^^ ■.k = 1,2,3} in voxel v are of primary interest and 
therefore treated as the treatment effects, whereas the differences of eigen- 
values between voxel v and its neighboring voxels {/Sj}"^]^ are the blocking 
effects and serve as nuisance parameters; {ej^k'-k = 1,2,3}"^]^ are random 
errors. Clearly, /3i = 0, such that (3.2) complies with (3.1) when j = 1. 

Nonetheless, as discussed above, based on the DT estimates, we can only 
obtain {Aj(fc),A: = 1,2,3}, the ordered version of {Xj^k^k = 1,2,3}. In other 
words, the values of Xj^k ^-re not available in practice. Therefore, the well- 
developed techniques for the additive model in a complete block design are 



LOCAL TESTS FOR ANISOTROPIC AREAS 9 

not directly applicable in analyzing model (3.2). We borrow its basic idea 
and integrate the corresponding contrast test statistic as one of the main 
parts in our proposed test statistic. The technical properties of our proposed 
test statistic in isotropically diffused DT voxels, that is, X%^ = \%-. = A^- 
then rely on the fact that Aj (fc) = A/^^n + f3j + £j^(k)- 

4. Local test based on neighboring eigenvalues. 

4.1. Test of hypotheses. Classifying the tissue types for a particular brain 
area plays a key role in the downstream fiber tracking algorithms. In voxels 
without fibers, water tends to diffuse in an isotropic manner, characterized 
by the eigenvalue property ALn ss X* ^ X*,. In contrast, water in voxels 
passed by fiber tracts usually presents anisotropic diffusion. DTs in these 
voxels typically have three possible morphologies, namely, prolate A/g-. > 
AL) ^ '^(i) 5 oblate AL^ ss ATg) > Ku ^i^d nondegenerate AL^ > ATj) > Ku j 
respectively, corresponding to brain voxels located in uniquely oriented fiber 
areas, crossed fiber areas with roughly the same fiber intensities and crossed 
fiber areas with distinct fiber intensities. In summary, anisotropy is char- 
acterized by either X*^^s > X%) or AL) > A?^n. Therefore, we consider the 
hypothesis testing problem, 

(4.1) /7o:A^3) = A^2) = A(i) vs. i7i:A^3)>A^2) or A^2) > A^^) 

for each voxel v. This testing problem can be reformulated as a form of 
contrast test, 

(4.2) Ho:AX* = vs. Hi:A{l,-)X*>0 for / = 1, or, . . . , or /i, 

where A* = {X%^, ^2)1^1))'^^ ^ is a full row rank matrix with rank(^) = fi 

and ^; ^i A[li,l2) = for /i = 1, . . . ,/i. There exists more than one possi- 
ble choice of A, such that hypothesis testing problems (4.1) and (4.2) are 
equivalent. For example. 



(4.3) A: 



1 -1 

1 -1 



which is adopted in our numerical studies. Clearly, /i = 2. 

To test (4.2), for each voxel v, we propose the test statistic K., represented 
by 

(4.4) IC = n(U„-0uJ^S-i(U„-0uJ, 



where U„ = (AA./^/MSE) • ^5.2; A. = (A.,(3), A.,(2), A.,(i))^; A.,(fe) = 
X]?=i \',{A:)/^5 A; = 1,2,3; S'^ is the mean of 5^ over the neighboring voxels 
of v; Sf= ELi_{h(k) - A„(.)} V2, i = 1, . . . , n; MSE = E"=i ZLiihik) " 
A.,(fc) — Aj (.) -|- A. (.)}2/{2(n — 1)}; let Vo be an estimate of Vq, the set of 
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isotropic voxels over the entire brain; 0u„ denotes the sample median of U„ 
over Vo; S denotes the sample covariance of y/nJJn over Vq- Vq is obtained 
from the following iteration steps: 

(i) Evaluate U.„'s over all voxels of interest. For some pre-given signifi- 
cant level a, let xfi-i-a be the (1 — a)th quantile of the chi-square distribution 
with n degrees of freedom. 

(ii) Let Vo,o include all voxels of interest. 

(iii) In the sth iteration (s = 1,2,...), let K. denote the statistic computed 
from (4.4) based on Vo,s-i, and let IK = cK. 

(iv) Include voxel v in Vo,s if the corresponding statistic K < Xh-i-o- 
(v) Repeat steps (iii) and (iv) until 0u„ converges. 

In step (iii) above, the constant c is to correct the bias in the evaluation 
of K. This is because in the sth iteration, Ou^ in IK is constructed by only 
using voxels in Vo,s-i, which includes voxels with IK capped by Xu,i-a ^^ the 
(s — l)th iteration. The value of c is derived as follows: 

E{Kl{ck < x%i.J} « -J-- Yl ^(U„ - ^uj"^S-i(U„ - ^uj 

I 1^0,; 



1 



|Vo,«-i| - 1 



■ Yl {^/^(u„-u.)fs-H^/^(u„-u.)} 

veVo,s-i 



|Vo,.-i| 

where "=" is followed by the fact that in step (iii), S is the sample covariance 
of y/nJJn. over Vo,s-i; I(-) is the indicator function. Therefore, 

(4.5) i?{IKI(IK < x%i-a)} = E{cKlicK < x^i-J} « c^^- 

On the other hand, according to Theorem 1 in Section 5, IK is approximately 
xfi distributed, 

(4.6) EmK<xl,.J]»l -j^-^^,"IUt, 

Combining (4.5) and (4.6), we set c = i /„""'" aw'rj./a) ^''^' *' 

Remark 1 . We now make some remarks concerning the construction of 
the statistic in (4.4). The main term U„ in IK consists of two parts: 



The first part y4A./vMSE mimics the t statistic of the contrast test for 
the additive model in a complete block design. Referring to the proofs of 
Theorems 1 and 2, under certain regularity conditions, it approaches infin- 
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ity with rate \/n, when there are significant differences among {A^^-, : k = 
1,2,3}, whereas it converges to a fixed constant when A?„n = AL\ = A^n. 
Therefore, it has good statistical power in identifying the differences among 
{A^^N ■.k = 1,2,3}, when model (3.2) is valid. 

• The second part y S'^ is added to increase the power of the test statistic 
K on the boundary of fibers. For any voxel on the boundary of fibers, 
in the sense that its selected neighborhood contains voxels belonging to 
both fiber and nonfiber areas, the assumption that the collected eigen- 
values in the selected neighboring voxels follow model (3.2) may not be 
appropriate. In this case, vMSE tends to in flate and, consequently, the 
statistical power of the first part AX./yMSE is limited. The second part 
then counteracts the effect of \/MSE, and therefore is particularly useful 
when the resolution of the DTI data is limited. 

• Furthermore, 6\j^, the sample median instead of sample mean of U^, is 
adopted just to ensure the robustness of the approach. 

Remark 2. In this paper we establish testing procedures for identifying 
the brain areas with fiber tracts. However, for voxels located in fiber tracts, 
their DTs have three possible morphologies, namely, prolate, oblate, and 
nondegenerate, respectively, corresponding to eigenvalue properties X%\ > 
AL^ ?» A^^ , A(3^ ^ X^-^ > A^^ , and A?3-j > X%^ > Xt. . One may also implement 
similar testing procedures as those in this paper to further classify these 
three possibilities. We leave the details out for presentational brevity. 

4.2. Adaptive selection of neighborhood. Following Sections 4.1 and 5, 
the asymptotic theories for K need that the neighborhood size n is large 
and that model (3.2) is satisfied. From the experimental point of view, as 
long as the resolution in a DWI experiment is sufficiently good, such that 
the proportion of voxels located on the boundary of fibers over the entire 
brain shrinks, we can simply employ a fix-shaped neighborhood in the con- 
struction of M. (e.g., choose the neighborhood as a fixed cube). However, for 
experiments with limited resolution, a fix-shaped neighborhood may not be 
a good choice, because in this case, the assumption that the eigenvalues in 
the neighboring voxels follow model (3.2) may not be well satisfied, in order 
to ensure n large required by Theorem 1. Such a problem is particularly 
severe for voxels located on the boundary of fibers. Therefore, development 
of a varied neighborhood is necessary. 

We propose an adaptive method to select the neighboring voxels based on 
the philosophy below. First, the adjacent voxels of v should have a better 
chance to be selected as neighboring voxels than those far away. Second, 
to ensure the validity of model (3.2), if the tensor in v is isotropic, the 
selected neighborhood should mainly consist of voxels with isotropic tensors. 
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Likewise, if the tensor in v is anisotropic, the selected neighborhood should 
be in favor of voxels with anisotropic tensors. Therefore, we incorporate the 
physical distances and similarity measures of DTs between voxel v and its 
nearby neighbors to establish the criteria for selecting neighboring voxels. 

For each voxel v and fixed number n of neighboring voxels, we summarize 
our proposed adaptive neighborhood selection approach as follows: 

(1) Fix a cube-shaped domain centered at v with reasonably large size 
X xy X z, whose voxels are candidates. Here x,y,z are integers and xyz > n. 

(2) Define the similarity score function / between voxel v and its neigh- 
boring candidate vi, 1 = 1,..., xyz, as 

f{v,vi) = dr,{'D{v),'D{vi))exp{C-dp{v,vi)}, 
where dp{v, vi) denotes the physical distance between voxels v and vi; dii{Y){v), 
■D(^O) — \/trace[{D(t;) — D(t>i)}^] is a measurement of the diffusion similar- 
ity between tensors in voxels v and vi [Alexander, Gee and Bajcsy (1999)]; 
C > is added to balance the contribution of dp{v,vi) and dD(D(w),D(u/)). 

(3) Select n voxels with the lowest / values as the neighboring voxels. 

We observe that C in / is adopted to balance the contribution of d£){'D{v), 
D(f;)) and dp{v,vi). When the resolution of the DTI data is high, dp{v,vi) 
is close to 0. Consequently, for any fixed C, exp{C ■ dp{v,vi)} approaches the 
constant 1. Therefore, for high resolution DTI data, the proposed approach 
is not sensitive to the choice of C. Throughout our numerical studies, we fix 
C = 0.1. 

We would like to point out that the proposed approach above is similar 
in spirit to the adaptive approaches in Tabelow et al. (2008) and Li et al. 
(2011), where iterative testing procedures are used to adaptively control the 
contribution of neighboring voxels in their proposed algorithms. Compared 
with their approaches, our approach is computationally more economic. The 
effectiveness of our proposed approach above has been demonstrated in Yu 
(2009) by simulation studies. 

We summarize our proposed procedure of constructing K in the supple- 
mental document [Yu et al. (2013)]. 

5. Theoretical properties. In this section we explore the theoretical prop- 
erties of our proposed test statistic K. The technical details are given in the 
supplemental document [Yu et al. (2013)]. Theorem 1 below establishes the 
asymptotic null distribution of K, when the number n of neighboring voxels 
is large. 

Theorem 1. Assume model (3.2) and Condition A in the supplemental 
document. Then for K. defined in (i-i); under the null hypothesis in (4-2), 
as n^ oo, 
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We would like to point out that the construction of K and the theoret- 
ical derivations of Theorem 1 are nontrivial and challenging. Following the 
discussion in Section 3, for each voxel v, we postulate that there exist un- 
observable one-to-one correspondences, that is, {(A^, A/^n) : A: = 1,2,3}, be- 
tween the estimated and true eigenvalues. The ordered eigenvalue estimates 
A(3) > A(2) > A(i), however, are neither unbiased estimates for X1^. > X1^. > 
Aqs, nor independent. We address these bias components in the construc- 
tion of K and the corresponding proof of Theorem 1 based on the intuition 
as follows. Referring to model (3.2), for an isotropic voxel v, the collected 
neighboring voxels can be modeled as Aj^(^.) = A^-, +/3j +ej (;,). As such, the 
bias components of eigenvalue estimates in the neighboring voxels of vjxre 
carried by ej^(fc), whose effects in our test statistic are counteracted by 0u„ 
constructed based on spatial information of the entire brain. 

To appreciate the discriminating power of K in the identification of aniso- 
tropic brain areas, the asymptotic power of K is established in Theorem 2 
below. 

Theorem 2. Assume model (3.2) and Condition A in the supplemental 
document. Then for voxel v, under the fixed alternative Hi in (4-2), as 
n — )• oo, 

where M is given by (A. 4. 3) in the supplemental document. 

Theorem 2 shows that as long as g(a(t;)) ^ g(b), Af > and IK — t- +oo 
at rate n, under the fixed alternative Hi. Here, g(-) is defined by (A. 3. 2) in 
the supplemental document; a.{v) = {E{2Si{v)}, E{Xi (3)(v)},£'{Ai (2)(^)}; 
S{Ai,(i)(^)})^;b = (£;(252),SK(3)(l)},SK(2)(l)},i?K(i)(l)}r;5f(^^) 
and 5^ are, respectively, the sample variances of {Ai (fc)(t;). A; = 1,2,3} and 
{^i (fc)(l),/c = 1,2,3}. Thus, under the fixed alternative, the power of our 
proposed test statistic K tends to 1 except in rare situations. Corollary 1 
below gives one specific example of M > 0. The proof is straightforward and 
omitted. 

Corollary 1. Assume conditions in Theorem 2. Suppose that Ej^ has 
a symmetric distribution about 0, that is, £j k has the same distribution as 
-Ej^k, and E{Xi^^^)} - E{Xi^(^2)} ^ E{Xi^^2)}'- E{Xi,^i)}. ThenM>0. 

6. Simulation study. 

6.1. Basic settings for numerical work. Since in a real DTI data set, the 
number of voxels, each of which corresponds to a hypothesis test, is typi- 
cally large, false discovery rate (FDR) techniques [Benjamini and Hochberg 
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Fig. 4. Left panel — neighbors of a voxel used in the FDR_t procedure. Right panel — gom- 
etry of the simulated brain. 

(1995), Storey (2002), Storey, Taylor and Siegmund (2004), Zhang, Fan and 
Yu (2011)] are incorporated in our numerical works to control the error 
rates. Two FDR procedures are employed in our study, namely, the conven- 
tional FDR procedure by Storey (2002) and the FDR^, procedure by Zhang, 
Fan and Yu (2011), which is capable of capturing the spatial information in 
imaging data. A short summary of the FDR/, procedure is provided in the 
supplemental document [Yu et al. (2013)]. 

Some settings of parameters throughout our numerical works are given as 
follows. For FDR and FDR^ procedures, set false discovery control level as 
0.01 and tuning parameter A = 0.2. The neighborhood for the FDR^ pro- 
cedure is set as the nearest 7 voxels shown in the left panel of Figure 4. 
When evaluating K for each voxel, A is given by (4.3). Apply the adaptive 
neighborhood selection approach proposed in Section 4.2 to select neighbor- 
ing voxels, where the domain for the candidate neighboring voxels is set as 
a 5 X 5 X 3 cube centered at v. There, n = 25 voxels are selected based upon 
/ with C = 0.1. The choice of n is referred to in the results in Section 6.3: 
when n = 25, the sampling distribution of K agrees reasonably well with the 
X^ distribution. 



6.2. Data simulation. We simulate several sets of DWI data over the 
entire brain. For each set, a 3D imaging space with the same brain ar- 
eas as the real DWI data in Section 7 is simulated, where fiber tracts are 
simulated to have the geometrical structure displayed in the right panel of 
Figure 4. DTs are simulated according to the locations of voxels in the imag- 
ing space. To this end, four distinct sets of eigenvalues, namely, [0.7,0.7,0.7], 
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[1.0,0.55,0.55], [0.8,0.8,0.5] and [0.9,0.7,0.5] (units: lO-^mmVs), are 
adopted to, respectively, simulate DTs with morphologies of isotropic (non- 
fiber areas), prolate (single blue and red bundles), oblate (intersected areas 
of blue bundles) and nondegenerate (intersected areas between blue and 
red bundles), such that all simulated DTs share the same mean diffusivity 
XJ-. = 0.7 X 10~'^mm^/s, a typical value in real human brains [Pierpaoli et al. 
(1996), Anderson (2001)]. 

Since the acquired attenuated signal intensity, (/)i(f), at each voxel v and 
gradient gj in real DWI data is typically generated by the square-root of the 
sum of squares of two random numbers in the DWI experiment [Henkelman 
(1985), Salvador et al. (2005), Zhu et al. (2007)], we simulate a reference 
signal {i = 0) and r = 12 diffusion attenuated signals (i = 1, . . . , 12) in each 
V [= {vx,Vy,Vz)] as 

Uv) = ^[(t>l{v)e^v{-bEj^*{v)^^i] + ei^^{v)f + ely{v), 

where (/>o(t;) = 1200 when v^ G (0,128], (l)l{v) = 1800 when v^, G (128,256]; 
the h factor h = 1000 when i > 0, 6 = when i = 0; diffusion gradients 
{gi : i = 1, . . . , 12} are adopted from the real DWI data in Section 7; D*(u) = 
Q{v)k*{v)Q'^{v)- A*(z;) = diag(A^3)(f),A*2)(t'),A*^)(7;)); Q{v) is a 3 x 3 or- 
thogonal matrix whose column vectors are composed of the eigenvectors of 
the simulated D*(f), 



Q{v) 



1/^/2 1/^/2 O' 

-1/^/2 1/^/2 

1 



or Q{v) 



1/^/2 -l/\/2 O' 

1/^/2 1/^/2 

1 



Clearly, the former corresponds to the red and the corresponding parallel 
narrow blue bundles in the right panel of Figure 4, while the latter models the 
other two blue bundles. The random errors £i^x{v) and Si,y{v) are simulated 
as independent and normally distributed with variance a^, which is varied to 
provide signal to noise ratios (SNRs), where SNR = (J)q[v)/(t. We examine 
four distinct SNRs, {5, 10, 15, 20}, each corresponding to one set of the 
simulated DWI data. 

6.3. Agreement between x^ distribution and K. With the DWI data sets 
simulated in Section 6.2, we estimate the DT in each voxel by regression 
model (2.1) and the corresponding eigenvalues by Schur decomposition. For 
each DWI data set, two sets (I and II) of K are constructed according to 
different settings of the adaptive neighborhood selection approach in Sec- 
tion 4.2. In particular, for Set I, the size of the candidate neighborhood 
and the number of selected neighboring voxels are, respectively, chosen as 
5x5x3 and n = 25, while those for Set II as 11 x 11 x 3 and n = 81. Other 
settings are given in Section 6.1. 
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n=25, SNR = 10 



n=25, SNR= 15 



n=25, SNR = 20 




Fig. 5. Empirical percentiles ofK (y-axis) versus percentiles of xfi distribution (x-axts). 
Top panels — candidate cubic 5 X 5 X 3, n = 25. Bottom panels — candidate cubic 11 x 11 x 3, 
n = 81. From left to right panels: SNR = 10, 15 and 20. Solid line — the 45 degree reference 
line. 

For each simulated data set, we collect all K's whose corresponding voxels 
are located inside the simulated nonfiber areas. The QQ plots of the (1st 
to 99th) percentiles of these K's against those of the xft distribution are 
displayed in Figure 5, with top panels based on Set I, bottom panels on 
Set II. The left, middle and right panels correspond to SNR = 10, 15 and 20, 
respectively. Results in Figure 5 demonstrate that the sampling distributions 
of IK, under both n = 25 and n = 81, agree reasonably well with the x^ 
distribution. 



6.4. Receiver operating characteristic curve. We compare our methods 
(K-FDR and K-FDRl) with FA with threshold (i.e., FA-threshold) and 
Smooth-FA with threshold (i.e. Smooth-FA-threshold) approaches by the 
receiver operating characteristic (ROC) curve, a widely adopted statistical 
tool for evaluating the accuracy of continuous diagnostic tests [Pepe (2003)]. 

Let {Ti, . . . , Tjv} be the set of statistics for all voxels over the entire brain, 
where Ty for any voxel v £ {!,..., A^} is the FA or Smooth-FA value, if 
the FA-threshold or Smooth-FA-threshold approach is used; is the p-value 
based on IK, if the K-FDR approach is used; is the p-value if the IK-FDR/, 
approach is used, where p stands for the median smoothed p-value of IK 
[Zhang, Fan and Yu (2011)]. For any given threshold t, if we classify a voxel 
V as anisotropic based on T^, G R{t), where R{t) = {Ty ■.Ty>t,v = l,..., N} 
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Fig. 6. ROC curves for SNR = 5, 10, 15 and 20. Dashed black line — p-value o/K; dotted 
blue line — p-value oflK; solid red line — FA; dash-dotted green line — Sniooth-FA. 



when r„ is the FA or Smooth-FA value; R{t) = {T„ -.T^, <t,v = 1, . . 
when Ty is the p- or p-value, then, 



.,A^} 



yN 



sensitivity: se(t) 



specificity: sp(t) 



EiUHTv(^R{t),Hiistrue} 



iVol 



^Af 



Z,UHTv^R{t),HoistTue} 

N-\Vo\ 



where |Vo| is the number of isotropic voxels over the entire brain. The ROC 
curve is then constructed as the 2D curve (se(t), 1 — sp(t)) when t ranges 
between and 1. The area under the ROC curve (AUC) is a popularly 
adopted measure of the accuracy of the test. More precisely, AUC is ranging 
between and 1, and the larger the AUC, the better the method. 

Following Sections 6.1 and 6.2, the ROC curves for FA, Smooth-FA, p 
and p are displayed in Figure 6 for data sets with SNR = 5, 10, 15 and 20, 
respectively. It has been seen from Figure 6 that the ROC curves of p (dot- 
ted blue lines) and p (dashed black lines) are consistently located above 
those of FA (solid red lines) and Smooth-FA (dash-dotted green lines) for 
all simulated data sets, indicating that K-FDR and K-FDR^ are capable 
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Table 1 
Sensitivity and specificity, SNR— 10, the FDR control level is 0.01 

Sensitivity Specificity 



FA > 0.3003 


0.8849 


0.4012 


Smooth-FA > 0.2803 


0.8852 


0.6242 


K-FDR 


0.7522 


0.9957 


K-FDRl 


0.8845 


0.9982 



of achieving better classification accuracy than FA-threshold and Smooth- 
FA-threshold approaches. Furthermore, the AUCs of p in all examined data 
sets are superior to those of p, suggesting that K-FDRj;, performs the best 
among all four approaches. 

6.5. Test results. In this section we present our test results in simulated 
data sets. Settings of our computations are given in Section 6.1. For the sake 
of clarity, we only present the results of SNR = 10. Those of SNR = 5, 15 and 
20 display similar phenomenon and are omitted. 

The FA threshold (>0.3003) and Smooth-FA threshold (>0.2803) are 
tentatively chosen such that their identified results share the same sensitivity 
as K-FDRi. The results of sensitivity and specificity by all four methods 
are displayed in Table 1. Clearly, K-FDR^ maximizes both the sensitivity 
(0.8845) and the specificity (0.9982) in this example. K-FDR achieves similar 
specificity (0.9957) but slightly smaller sensitivity (0.7522) compared to K- 
FDR^. However, in order to yield comparable sensitivity (0.8845) as K- 
FDRi, FA-threshold and Smooth-FA-threshold approaches produce much 
lower specificities (0.4012 for FA; 0.6242 for Smooth-FA). 

The performances of the four methods are further compared on two se- 
lected axial slices. Throughout our simulation and real data examples, we 
apply the same registration transformations from the brain data to the Tl 
high-resolution image of the subject's brain. The two slices with the simu- 
lated brain anisotropic areas highlighted are given in the leftmost panel of 
Figure 8. Figure 7 displays the color maps of the FA, Smooth-FA, — log(p) 
and — log(p), with all — log(p) and — log(p) values greater than 10 set 
equal to 10 to improve the visualization. Figure 8 compares the detected 
anisotropic areas by FA > 0.3003, Smooth-FA > 0.2803, IC-FDR and K- 
FDRi for SNR = 10. 

As clearly evidenced in Figures 7 and 8, K-FDR and K-FDR/, not only 
provide detected results with better accuracy than FA-threshold and Smooth- 
FA-threshold approaches, but also yield better contrasts between the signif- 
icant areas and the nonsignificant ones. In contrast, the results from FA 
(>0.3003) and Smooth-FA (>0.2803) not only fail to detect some truly sig- 
nificant voxels, but also present tiny scattered faulty findings, which expect 
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Fig. 7. Comparison of color maps for simulated data set. From, left to right: FA; 
Smooth-FA; - log(p) by K; -log(p) by K. SNR= 10. 

to contaminate the downstream fiber tracking results. It has been seen in 
the second to right and rightmost panels of Figures 7 and 8 that the de- 
tected results by K-FDR and IK-FDR^ well capture the primary features 
of the simulated anisotropic areas. Compared with IK-FDR, the IK-FDRj;, 
approach offers slightly more accurate identifications in both isotropic and 
anisotropic water diffusion areas. 




Fig. 8. Comparison of bram anisotropic areas discovered for the simulated data set. 
From left to right: simulated brain anisotropic areas; FA > 0.3003; Smooth-FA > 0.2803; 
K-FDR; K-FDRz,. SNR= 10. The control level is 0.01. 
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7. Real data example. We apply our proposed testing procedures on 
five subjects, whose DWls were acquired by the magnetic resonance (MR) 
experiments described below. 

The brain magnetic resonance images (MRIs, including DWI and iMRI) 
of each subject were acquired with a GE SIGNA 3-T scanner equipped with 
high-speed gradients and a whole-head transmit-receive quadrature birdcage 
headcoil (GE Medical Systems). The anatomical scan for each subject took 
approximately 20 minutes [Dalton et al. (2005)]. In the anatomical scan- 
ning, the size of each voxel in an xy-plane is 0.9375 mm x 0.9375 mm, field 
of view = 24 cm^, matrix = 256 x 256; 30 axial slices are acquired along the 
z-axis, slice thickness = 3 mm. A single reference image at 6 = and 12 
diffusion-attenuated images with noncoUinear directions of diffusion gradi- 
ents at 5 = 1000 s/mm^ were obtained. Since we focus on the analysis of 
the anatomical structures of the human brain in this paper, the detailed 
information for the functional scans is omitted. 

Using the DWI data of a single subject as the representative, we first 
present and compare the results by all four methods on two selected axial 
slices of the brain. The results for the other four subjects display similar 
scenarios and are omitted. 

The acquired data set contains 256 x 256 x 30 = 1,966,080 voxels with 
400,309 voxels located inside the brain. In each voxel, the DT is estimated 
from regression model (2.1). After that, the corresponding eigenvalues are 
obtained by Schur decomposition. 

All settings are given in Section 6.1. The color maps of FA, Smooth- FA, 
— log{p) and — log(p) are displayed in Figure 9 on two selected axial slices. 



V . 






■1 


m 


m: 






■V 




If 


D 1 Ci -: 





Fig. 9. Comparison of color maps for the real bram data set. From, left to right: FA; 
Smooth-FA; -log(p) hyK; -log(p) hyK. 
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Fig. 10. Comparison of brain anisotropic areas discovered for the real brain data set. 
From left to right: FA > 0.35; Smooth-FA > 0.35; -log(p) by K; -log(^ by K. The 
control level is 0.01. 



whereas the corresponding detected anisotropic diffusion brain areas by ah 
four methods are provided in Figure 10. As evidenced in Figure 10, compared 
with the identified anisotropic areas by K-FDR or IK-FDR^, FA-threshold 
and Smooth- FA-threshold approaches produce more noisy detections. For 
example, inspection of areas highlighted by red rectangles in the top panels 
of Figure 10 (enlarged in Figure 11), FA > 0.35 and Smooth-FA > 0.35 detect 
more scattered tiny areas than K-FDR and K-FDR^. Those are highly likely 
to be faulty findings. Furthermore, an overview of the areas located close to 
the top of the highlighted areas, FA > 0.35 and Smooth-FA > 0.35 present 
more findings than K-FDR and IK-FDR^. However, those areas are located 
close to the boundary of the brain, and are most likely to be nonfiber areas. 
In the meantime, as illustrated in Section 2.3, since the number of gradients 
r in each voxel is small, FA-threshold and Smooth-FA-threshold approaches 
cannot clearly infer and control the error rate of the identification, and 




Fig. 11. Enlarged areas in the rectangles of Figure 10. Rotated 
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Table 2 
Number of voxels earned in Si,m (le.ft) and S2,m (right) over each subject's brain 









Subject 










Subject 






Methods 


1 


2 


3 


4 


5 


1 


2 


3 


4 


5 


FA > 0.35 


323 


360 


348 


356 


280 


718 


770 


720 


758 


671 


Smooth-FA > 0.35 


279 


288 


327 


298 


227 


775 


793 


657 


721 


660 


K-FDR 


155 


177 


166 


173 


189 


315 


409 


357 


403 


355 


K-FDRi 


63 


77 


100 


93 


103 


129 


158 


149 


212 


138 



therefore lack the rigorous criterion in selecting the appropriate thresholds 
in practice. The superiority of our proposed methods to FA-threshold and 
Smooth-FA-threshold approaches can be further illustrated by Figure 9, the 
color maps of FA, Smooth-FA, — log(p) and — log(p). The color maps of 
both — log(p) and —log{p) show better contrasts between the anisotropic 
and isotropic areas than those of FA and Smooth-FA, indicating that K.- 
FDR and IK-FDR^ more effectively separate anisotropic diffusion areas from 
isotropic ones. 

To further illustrate the efficacy of our methods in reducing the scattered 
faulty findings, we summarize the frequency of identified isolated anisotropic 
voxels over each subject's brain for all four methods in Table 2. In particu- 
lar, denote by M{v) = {v' = {v'^, Vy, v'^) :\v[ — Vi\<l, for all i = x,y, and z} 
the nearest neighboring voxels of v. Table 2 displays the number of voxels 
carried in Si^m (left) and 52,m (right) over each subject's brain. Here Su,m = 
{v. method m identify u voxels mAf{v) as anisotropic}. Clearly, Si^m car- 
ries voxels where only the voxel v itself is identified by method m over voxels 
in J\f{v). Likewise, 52, m contains voxels where only two voxels, that is, the 
voxel V itself and another voxel, are identified by method m over voxels 
in M{v). We observe that voxels in Si^m and 52,m are highly likely to be 
faulty findings by method m, since they are "isolated" from other identified 
anisotropic voxels; fiber tracts, in contrast, are typically spatially connected. 

It has been seen from Table 2 that our methods continue to outperform 
FA-threshold and Smooth-FA-threshold approaches by producing much less 
isolated findings. Compared with K-FDR, IK-FDR/, produces even a smaller 
number of isolated identifications. Such a result is not surprising considering 
what has been observed in Figures 7-10. Combining all the numerical results 
above, we therefore recommend the identifications by K-FDRj;, as the final 
results. 

8. Discussion. In DTI studies, one of the important research topics is 
to refine the identification of the anisotropic water diffusion areas of human 
brain in vivo. There are two general strategies aiming to address this prob- 
lem. The first is to improve the DT estimation. A downstream procedure is 
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then needed to identify the anisotropic water diffusion areas. The second is 
to refine the construction of the scalar measurements or estabhsh more pow- 
erful test statistics for every brain voxel. The identification is then based on 
thresholding the measurements or a certain testing procedure. We observe 
that the second provides more intrinsic insight into the water diffusivity in 
each voxel, and therefore is more effective in allusion to the identification of 
anisotropic water diffusion areas. 

From an experimental point of view, there are two ways to improve the 
acquisition schemes for the DWI data. One is to increase the number of diffu- 
sion gradients for every brain voxel, the other is to improve the resolution of 
the imaging space, that is, increase the number of brain voxels. To the best 
of our knowledge, existing methods for constructing scalar measurements 
or test statistics are all single voxel based. Therefore, the corresponding in- 
ferences improve only when the number of diffusion gradients in each voxel 
increases, while ignoring the possible improvement of the resolution of the 
imaging space. The methods proposed in this paper fill this gap by incorpo- 
rating the eigenvalues in the neighboring voxels in the construction of the 
test statistic K. 

In this study we have established the asymptotic distribution of our pro- 
posed test statistic. One of the main assumptions required by our theoretical 
results is that the number of neighboring voxels for constructing IK is large. 
This assumption can be well achieved when the resolution of the imaging 
data is high. As such, the bias components carried in the eigenvalue estimates 
no longer play a key role in the identification of anisotropic water diffusion 
brain areas. In both simulation and real data analysis, we have observed 
that our proposed EC-FDR and K-FDRj;, approaches lead to different iden- 
tification results from FA-threshold and Smooth-FA-threshold approaches, 
popularly adopted in the DTI community. In particular, the scattered find- 
ings by our methods are much less than those by FA-threshold and Smooth- 
FA-threshold approaches, indicating that by incorporating neighboring in- 
formation, our methods are capable of screening out those isolated voxels 
which are highly likely to be faulty findings. Results based on simulated 
DWI data demonstrate that our proposed test statistic IK agrees reasonably 
well with the x^ distribution when n = 25 (or larger), and our methods 
achieve better accuracy than FA-threshold and Smooth-FA-threshold ap- 
proaches in the identification of anisotropic brain voxels. Furthermore, the 
Smooth-FA-threshold approach is capable of partially solving the bias prob- 
lem in the eigenvalue estimates [Polzehl and Tabelow (2009)]. However, the 
performance of the approach heavily replies on the estimation of the het- 
eroscedastic variances over the entire brain. These variances, in turn, are 
modeled by a linear model and estimated using the reference signals (/>o('v). 
We observe based on simulation studies that when the reference signals over 
the entire brain are not homogeneous or do not share comparable variances 
as attenuated signals (j)i{v), the performance of the approach varies [see 
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more simulation results provided in the supplemental document: Yu et al. 
(2013)]. In contrast, under all these cases, our proposed K-FDR^ approach 
consistently offers descent results. We therefore conclude that over all four 
methods, our proposed IK-FDR^ approach performs the best over all our 
simulation studies. Unlike the simulation examples, we are unable to show 
the true anisotropic (fiber) areas of the human brains in real DTI data. 

We would like to point out that in DTI studies, identification of anisotropic 
water diffusion areas is just one step of the full analysis. Downstream anal- 
ysis, such as fiber tracking, is usually needed to fully capture the physical 
structure of the human brain. 

In this paper we have focused on establishing testing procedures to distin- 
guish anisotropic DT voxels from isotropic ones based on second order DT 
models. Our proposed methods have been compared with the FA-threshold 
and Smooth-FA-threshold approaches. The results of this paper can serve as 
a benchmark in analyzing the anatomical structures of human brains based 
on DTI data. We observe that the following topics are highly related to this 
paper, and may absorb the interest of the community in future research: 

• With the idea of incorporating spatial information, single-voxel-based ap- 
proaches, including FA, can possibly be improved by appropriately ac- 
counting for the information in the neighboring voxels. Furthermore, in 
this paper, K is established on the DTs estimated from model (2.1). Sim- 
ilar approaches can be constructed based on more sophisticated DT esti- 
mates from other approaches. 

• A similar strategy as that in this paper can be adopted to establish test- 
ing procedures for teasing apart the morphologies of anisotropic DTs. 
Furthermore, the basic ideas can be further extended to establish testing 
procedures for identifying the presence of signals for data with spatial 
structures, though we focus on DTI data in this paper. 

• Another popular topic in DTI research is to consider higher-order tensor 
models [Grigis et al. (2011) and therein], which are powerful tools for 
investigating the fiber structures when there are several fiber bundles in 
a single voxel. The local test idea in this paper can be possibly extended 
to identify the number of intersected fiber bundles based on those higher- 
order tensor models. 

Acknowledgments. We thank the Editor, the Associate Editor and four 
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SUPPLEMENTARY MATERIAL 

Supplement to "Local tests for identifying anisotropic diffusion areas in 
human brain with DTI" (DOI: 10.1214/12-AOAS573SUPP; .pdf). This file 
provides proofs for Theorems 1 and 2, a short summary of FDR^ procedure 
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[Zhang, Fan and Yu (2011)], steps for constructing K based on DWI data 
and some more simulation results. 
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